Why use Mizer?

Mizer is a tool that can be used to simulate a dynamic size spectrum in an marine ecosystem, subject to changes through time, such as fishing pressure. Multiple interacting species and fishing gears can be defined allowing for a range of fisheries management strategy scenarios, and their ecosystem impacts, to be tested.

Below we explore some examples of what can be done using a previously calibrated Mizer model.

Species Collapse and Recovery in a Heavily Fished Multispecies System

In this example we will read in a calibrated model for the North Sea, which consists of 12 interacting species and a background resource. All species are predators and prey since they feed according to size-based rules and they encounter each other. Fishing is applied as a logistic curve by length, matched to the maturity length \(l_{50}\) of each species.

More information on how this model was calibrated with historical data can be found in our “mizerHowTo” package (tutorial("HTM3")) and in this blogpost.

First let’s take a look at what the model has used in the historical period (from ICES stock assessments).

The species’ maximum fishing mortality rates (which we call for convenience effort) for Cod and Saithe have declined but are not as low as what would be in line with single-species “Fmsy” (a reference level of fishing deemed sustainable from single-species stock assessements) and the fishing mortality rate for Sprat and Dab is very high.

Then you can take a look at how the modelled species biomasses change through time by running the following:

plotlyBiomass(sim, start_time = 1950, end_time = 2018)
# size spectra averaged for the final 5 years
plotlySpectra(sim,time_range = c(2015:2018),power=2)

Examining changes in the community relative to an unfished state

To be able to examine any changes in the community relative to an unfished state we can re-set effort = 0 and calculate the steady state.

sim0 <- projectToSteady(params, effort = 0, return_sim = TRUE, t_max = 200)
plotlyBiomass(sim0)

We can compare the current size spectra (averaged over 2015-2019) to the unfished size spectra to assess whether there is any evidence of a size-structured trophic cascade due to fishing.

plot_relative_biomass(sim0, sim)

Here we can see the effect of the reduction in large sized individuals of heavily fished species on the other sizes and species in the model, relative to the unfished steady state. The abundance of some (but not all) of the smaller to medium sizes of prey are a lot higher when their larger predators are removed (note the logarithmic scale). Sprat looks to be much lower.

We can do the same with Biomass to see when, if any, of the species collapse. For simplicity, we use < 0.1 of B/B_unfished as a proxy for a reference point for population collapse. We can add other reference points to this kind of plot, for example a simple rule of thumb for B_msy could be 0.5*B_unfished.

Projecting future recovery scenarios

We may now wish to explore the potential recovery of the larger species and sizes in the system. To do this we set up another scenario, where the model starts with the last time step of the fished scenario.

The effort values in 2018 were

Current effort
Sprat Sandeel N.pout Dab Herring Gurnard Sole Whiting Plaice Haddock Saithe Cod
1.24 0.0845412 0.311 1.278 0.191 0.3773742 0.335 0.198 0.192 0.247 0.419 0.3200014

Let’s start a new simulation that begins with the effort from 2018 and projects forward for 50 years. We will apply a linear reduction in effort for a selected species to a target value (here assumed for simplicity to be 0.2 with effort expressed in terms of the species’ fishing mortality rate for fully selected sizes).

To do this we need to work with the effort array (time x gear) to enable changes in effort through time, for this scenario. Here we have a ‘gear’ for each species, since effort in this case were annual single-species fishing mortality estimates.

We can see that when we reduce fishing on Sprat it increases the biomass of this species but also affects the biomass of other species in the community.

Are any species collapsed still?

Example 2: Set up your own fishing scenario

Rather than an entire time-series, we can also simply examine differences between two time-averaged states under different fishing regimes.

We can alter the fishing parameters using a function called gear_params() and by changing the effort input.

Let’s take a look at the fishing parameters. Note that the gears in the above model were already setup to be species-specific.

gear parameters
species gear sel_func l25 l50 catchability initial_effort
Sprat Sprat sigmoid_length 7.6 8.1 1 1.2953333
Sandeel Sandeel sigmoid_length 9.8 11.8 1 0.0651055
N.pout N.pout sigmoid_length 8.7 12.2 1 0.3138000
Dab Dab sigmoid_length 11.5 17.0 1 0.9780000
Herring Herring sigmoid_length 10.1 20.8 1 0.1815000
Gurnard Gurnard sigmoid_length 19.8 29.0 1 0.4625057
Sole Sole sigmoid_length 16.4 25.8 1 0.3738333
Whiting Whiting sigmoid_length 19.8 29.0 1 0.2426667
Plaice Plaice sigmoid_length 11.5 17.0 1 0.1848333
Haddock Haddock sigmoid_length 19.1 24.3 1 0.3015000
Saithe Saithe sigmoid_length 35.3 43.6 1 0.3930000
Cod Cod sigmoid_length 13.2 22.9 1 0.2666675

We can group these species together according to the gears they are caught by. Let’s imagine a big super trawler.

# allocate species to gear types
gear_params(params)$gear <- c("super_trawler")

Note that catchability is set to 1. This is because the fishing “effort” was here assumed to be the fishing mortality rate of fully selected sizes (see here setFishing).

The previous effort won’t work with these new gears, as it is gear x time. We only have a single gear now, so this is easier to set up.

params <- setFishing(params, initial_effort = 1)

Now let’s run two simulations, one with light fishing effort (effort = 0.5) and one heavely fished (effort = 1.5).

sim_light <- project(params, effort = .5, t_max = 100)
sim_heavy <- project(params, effort = 1.5, t_max = 100)
plot_relative_biomass(sim_light,sim_heavy)

Large individuals are the most affected by the high fishing effort. This create a trophic cascade where smaller individuals, even if they are affected by the same fishing effort, will increase in abundance due to the release in predation from larger individuals.

Exercise: Now try to edit the above code to create your own fishing scenario.

Example 3: add a species to the North Sea model

Let’s start with the North Sea model as in Example 1.

params@species_params$biomass_observed = getBiomass(sim)[nrow(sim@n), ] # set biomass_observed

plotCalibration(sim)

plotGrowthCurves(sim, species_panel = TRUE)

That is set up and behaving OK. Now let’s add another species and see what happens. We call it ‘Colossal squid’, arbitrarily pick the species parameters, and an observed biomass of 110g.

species_params <- data.frame(
  species = "Colossal squid",
  w_inf = 3e4,
  w_mat = 2e3,
  beta = 100,
  sigma = 2,
  k_vb = 0.2,
  a = 0.01,
  b = 3
)

new_params = addSpecies(params, species_params)
## No h provided for some species, so using f0 and k_vb to calculate it.
new_params@species_params$biomass_observed[13] = 1e10 # set observed biomass
new_params = steady(new_params)
## Convergence was achieved in 67.5 years.
sim2 = project(new_params, t_max = 5)
plotlyBiomass(sim2, power = 2)
plotBiomassObservedVsModel(sim2)

Biomass is way off for squid! Use Gustav’s method in this blog post for retuning parameters.

new_params <- new_params |> calibrateBiomass() |> matchBiomasses() |> steady() |>
  calibrateBiomass() |> matchBiomasses() |> steady() |>
  calibrateBiomass() |> matchBiomasses() |> steady()
## Convergence was achieved in 73.5 years.
## Warning in setBevertonHolt(params): For the following species `erepro` has been
## increased to the smallest possible value: erepro[Colossal squid] = 1.97e-07
## Convergence was achieved in 30 years.
## Warning in setBevertonHolt(params): For the following species `erepro` has been
## increased to the smallest possible value: erepro[Colossal squid] = 1.47e-07
## Convergence was achieved in 4.5 years.
plotBiomassObservedVsModel(new_params)

Done now! Look at the steady state size spectrum.

plotlySpectra(new_params, power = 2)

And look at the growth curves, compared to von Bertalanffy growth curves.

plotGrowthCurves(new_params, species_panel = TRUE)

Have a peek at the reproductive parameters, to check.

species biomass_observed R_max erepro
Sprat Sprat 9.703980e+10 2.038528e+10 0.0097946
Sandeel Sandeel 4.636213e+12 7.710521e+13 0.0000719
N.pout N.pout 3.005304e+11 2.371267e+10 0.0070738
Dab Dab 6.764875e+11 2.032365e+10 0.0105072
Herring Herring 8.989985e+10 1.985117e+09 0.0114390
Gurnard Gurnard 7.067464e+10 1.833829e+09 0.0002563
Sole Sole 3.572392e+10 3.360454e+08 0.0098435
Whiting Whiting 9.066888e+10 1.132840e+09 0.0033342
Plaice Plaice 2.397591e+11 8.976962e+09 0.0017399
Haddock Haddock 2.378042e+11 1.304191e+09 0.0029266
Saithe Saithe 1.318700e+11 1.264573e+08 0.0636085
Cod Cod 5.515711e+11 7.844713e+08 0.4499833
Colossal squid Colossal squid 1.000000e+10 Inf 0.0000001